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Background: The linear response of the nucleus to an external field contains unique information about the 
effective interaction, correlations governing the behavior of the many-body system, and properties of its excited 
states. To characterize the response, it is useful to use its energy-weighted moments, or sum rules. By comparing 
computed sum rules with experimental values, the information content of the response can be utilized in the 
optimization process of the nuclear Hamiltonian or nuclear energy density functional (EDF). But the additional 
information comes at a price: compared to the ground state, computation of excited states is more demanding. 

Purpose: To establish an efficient framework to compute energy-weighted sum rules of the response that is 
adaptable to the optimization of the nuclear EDF and large-scale surveys of collective strength, we have developed 
a new technique within the complex-energy finite-amplitude method (FAM) based on the quasiparticle random- 
phase approximation (QRPA). 

Methods: To compute sum rules, we carry out contour integration of the response function in the complex- 
energy plane. We benchmark our results against the conventional matrix formulation of the QRPA theory, the 
Thouless theorem for the energy-weighted sum rule, and the dielectric theorem for the inverse energy-weighted 
sum rule. 

Results: We derive the sum-rule expressions from the contour integration of the complex-energy FAM. We 
demonstrate that calculated sum-rule values agree with those obtained from the matrix formulation of the QRPA. 
We also discuss the applicability of both the Thouless theorem about the energy-weighted sum rule and the 
dielectric theorem for the inverse energy-weighted sum rule to nuclear density functional theory in cases when the 
EDF is not based on a Hamiltonian. 

Conclusions: The proposed sum-rule technique based on the complex-energy FAM is a tool of choice when 
optimizing effective interactions or energy functionals. The method is very efficient and well-adaptable to parallel 
computing. The FAM formulation is especially useful when standard theorems based on commutation relations 
involving the nuclear Hamiltonian and external field cannot be used. 

PACS numbers: 21.60.Jz, 21.10.Re, 23.20.Js, 24.30.Cz 


I. INTRODUCTION 

Atomic nuclei exhibit various kinds of collective exci¬ 
tations, with characteristics considerably different from 
simple nucleonic excitations mm- Among those, giant 
resonances form a distinct class [3]. Although their exci¬ 
tation energies are relatively high compared to the low- 
energy collective modes, the main characteristics of giant 
resonances are understood in terms of the superposition 
of many nucleonic excitations. Experimentally, various 
types of giant resonances have been seen. Examples are 
shape vibrations, spin excitations, and charge-exchange 
excitations of various multipolarity and isospin. These 
modes carry rich information about basic nuclear prop¬ 
erties. 

There has been excellent progress in the modeling of 
atomic nuclei using nuclear density functional theory 
(DFT) g], State-of-the-art energy density functionals 


(EDFs), optimized to various classes of data gHH] en¬ 
able a quantitative description of global nuclear proper¬ 
ties throughout the nuclear landscape [TOUTS] . Ground- 
state properties of nuclei, such as binding energies, charge 
radii, effective single-particle energies of doubly-closed 
shell nuclei, and basic parameters characterizing the nu¬ 
clear matter equation of state, are typically used as em¬ 
pirical inputs in EDF parameter optimization. However, 
properties of excited states, such as giant resonances, are 
seldom considered (see Refs. Eonra for representative 
examples of work along those lines). This results in large 
uncertainties of EDF parameters sensitive to, and gov¬ 
erning, low- and high-frequency nuclear excitations. The 
EDFs of the next generation are expected to overcome 
this deficiency by including selected properties of the gi¬ 
ant resonances into the pool of observables used in the 
optimization. 

To extract the information content of giant resonances, 
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the sum rule technique has been widely used. 

For instance, mean giant resonance energies can be re¬ 
lated to the ratio of the sum rules of different energy 
moments [221 125H28] • The inverse energy-weighted sum 
rule provides information about the nuclear polarizabil¬ 
ity, which is the fundamental quantity characterizing the 
nuclear response. An important quantity, in the context 
of studies of neutron-rich matter, is the electric dipole 
polarizability, which is related to symmetry energy and 
its density-dependence j29H3T] . Various polarizabilities 
carry information about instabilities in nuclear matter 
[32H35. In some cases, the Thouless theorem [5oH5T] 
provides a simple way to access sum rules directly from 
the Hartree-Fock-Bogoliubov (HFB) solution. Unfortu¬ 
nately, the Thouless theorem applies to positive-odd en¬ 
ergy moments, and simple expressions can be derived 
only for simple operators (such as multipole moments). 
Moreover, the theorem is justified only if a Hamilto¬ 
nian representation of the interaction is available, which 
is generally not the case for nuclear DFT where mod¬ 
ern EDFs are usually not connected to an underlying 
Hamiltonian, and often break local gauge invariance [58] . 
Therefore, an efficient technique to compute nuclear sum 
rules, regardless of the form of the operator F, is desired. 

The direct evaluation of sum rules from self-consistent 
QRPA matrix solutions is computationally demanding 
because of configuration spaces involved. A recent for¬ 
mulation of the sum rule in terms of QRPA matrices 
enables the computation of sum rules without diagonal¬ 
izing the QRPA matrix m- Nevertheless, this method 
still requires knowledge of the QRPA matrix, which has 
large dimensions, especially when spherical symmetry is 
broken. Other recent developments include applications 
of the Lanczos algorithm to RPA sum rules [55] and the 
use of the Lorentz integral transform method and the 
Lanczos technique [40] . 

The finite amplitude method [41], based on the linear- 
response approach, significantly reduces the computa¬ 
tional cost of the QRPA problem. The residual two- 
body interaction is numerically computed from the finite- 
amplitude nucleonic fields induced by an external po¬ 
larizing field. The FAM has been recently imple¬ 
mented to various self-consistent frameworks, including 
three-dimensional HF [41], spherical HFB [42], axially- 
deformed Skyrme-HFB |45M5j . and relativistic mean- 
field models [HI 07]. The FAM has been applied to 
the description of giant resonances and low-energy dipole 
strength [48] |49] , the computation of the QRPA matrix 
elements [50], and the description of discrete low-lying 
QRPA modes by means of the contour integration tech¬ 
nique in the complex energy plane m- 

The objective of this study is to propose an efficient 
approach to sum rules by using the contour integration 
technique of Ref. ED- Because of its inherently paral¬ 
lel structure, the new method is ideally suited for op¬ 
timizations of next-generation nuclear EDFs, informed 
by experimental data on multipole and charge-exchange 
strength. This paper is organized as follows. Section [IT] 


summarizes the basic expressions. In Sec. |III[ we present 
the formulation of the complex-energy FAM approach to 
sum rules. Section |IV| contains numerical tests, bench¬ 
marking examples, and applications to realistic cases. 
The conclusions and outlook are given in Sec. [V] 


II. BASIC EXPRESSIONS 
A. Sum rule 

The ground-state (g.s.) strength function S(E) for a 
one-body operator F is defined as 

S(E) = Y / S(E-E v )\(u\F\0)\ 2 , (1) 

V 

where |0) and \v) denote, respectively, the ground state 
and excited state of the system. The fc-th moment of 
S{E), 


m k (F) = J(E- E 0 ) k S(E ) dE, (2) 

is called the energy-weighted sum rule of order k. In 
terms of the transition matrix elements of F, it is given 
by: 


m k {F) = Y,{E v -E 0 ) k \{v\F\Q>)\\ (3) 

V 


As discussed in, e.g., Refs. mm, certain sum rules are 
independent of the specific many-body theory used to 
describe the ground state and excited states. For exam¬ 
ple, the nuclear shell model and QRPA frameworks have 
been widely used to evaluate the sum rules. In QRPA, 
the excitation energy E u — Eq is replaced with the QRPA 
frequency fl„, which is the eigenvalue of the matrix equa¬ 
tion: 



(4) 


where A and B are QRPA matrices. The QRPA equation 
Q has positive-energy solutions VL„ > 0 (v > 0) with 
(X l ',Y v ), and mirror negative-energy solutions = 
-f x < 0 with {X~ v , Y~ v ) = ( Y v *,x v *). The positive 
frequency solutions, being the physically relevant ones, 
are used for the sum rule, and the summation in Eq. |3]) 
is, therefore, restricted to QRPA modes with v > 0. 


B. Finite amplitude method 

The FAM is an efficient technique to obtain the re¬ 
sponse function S(E) without explicitly computing the 
A and B QRPA matrices in Eq. For the details per¬ 
taining to FAM, we refer the reader to, e.g., Ref. 02] , 
The complex response function for a given operator F 
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at a given complex frequency w 7 = ui + iq, found as a 
solution of the FAM equations, is given as 
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E 

i />0 


l(E|o)l 2 | l(Q|Ak)l 2 

fly CU 7 fly I CU 7 


(5) 


The Lorentzian distribution of the strength function is 
obtained by taking the imaginary part of S: 

-Im S(F, w 7 ) 

7T 

_7 V [ kEiq)i 2 kqiE)i 2 1 

7T^ 0 \(^-a;)2+ 7 2 (n v + w ) 2 + 7 2 j • U 

A contour integration along the path C v , which encircles 
a real energy pole fly of the response function, gives the 
QRPA transition strength to state \v) m- 


where the contour D encircles all the positive QRPA fre¬ 
quencies fly > 0, and excludes all the singularities of the 
complex function f(u>). By setting f(u>) = w fc , we obtain 
the expressions for the sum rule rrik(F). 

In the following, we assume the operator F to be Her- 
mitian for simplicity. In this case, positive and negative 
energy solutions are associated with the same transition 
strength: 

IHF|o)| 2 = |(o|i>)| 2 . (ii) 


The above equation does not hold when F is not Hermi- 
tian. However, Eq. (10) still can be used with an appro¬ 
priately chosen contour D. 


A. Laurent series of the FAM response function 


E/ S(F,^)ckj 7 = \(u\F\0)\ 2 (fly>0), (7) 

J Cv 

or, alternatively, along C_y, 

E£ S(F,u>^)duj^ = — |(0|F|^)| 2 

= -|< E t | 0)| 2 ( fi -«,< 0 ). ( 8 ) 

For a small 7 <C w, the relation l/(w + * 7 ) = P(1/oj) — 
inS(u}) holds, and the sum rules can be formally calcu¬ 
lated using 

A 1 r°° 

rrik(F) = -lim / w fe Im S(F, u) + i^)duj . (9) 

7T 7->0 J 0 

An approximate value of the sum rules can be found from 
this expression from a finite value of 7 pl2] 1431 2ES] . 
However, to guarantee sufficient numerical accuracy, a 
very fine mesh would be required for the integration to 
take into account all the QRPA modes, whose locations 
are not known beforehand. 


III. SUM RULE EXPRESSIONS IN FAM 

In this section we introduce the sum rule approach 
based on the contour integration of the FAM. For sim¬ 
plicity, we assume that the operator F cannot excite spu¬ 
rious modes, and all the QRPA energies fly are non-zero. 
We also assume that the HFB state is stable with respect 
to small density variations, i.e., there are no imaginary- 
frequency QRPA solutions. This guarantees that all the 
QRPA poles lie on the real axis. In the following, we 
shall adopt the notation ui for a complex frequency. 

The basic idea behind the FAM approach to sum rules 
is to utilize the identity based on Cauchy’s integral the¬ 
orem: 

f(u } )S(F,co)dco = V/(LQ)|HF|0)| 2 , (10) 

Jd „>o 


By using the Laurent series expansion of (1 — z)~ x , we 
can derive the expansion of the FAM response function. 
The FAM response function has poles at u> = and 
—fIn the inner region below the lowest QRPA pole, 
|<j| < min ;/> o fli/, S(F,u>) can be written as 


S(F,u } ) = ~2j2m- {2 n+i)( P )“ 2n - ( 12 ) 

n—0 

One can see that odd-fc sum rules can be simply related 
to the expansion coefficients of The same is true 

in the outer region above the highest QRPA pole, oj > 
max„>o fly, where the response function can be expanded 
as 


s(P,.)=2f:^Fn. 


(13) 


n —0 


The expansions (12 1 and (131 are generalizations of 
expansions proposed in Ref. m to the full complex en¬ 
ergy plane. We note that the inverse energy-weighted 
sum rule (k = — 1 ) is found by setting ui = 0 in 
Eq. (12). This should be done with care, however. If 


spurious modes are present, they would produce a zero- 
frequency pole resulting in numerical instabilities near 
or at the pole. If we choose the semicircle Ai (counter¬ 
clockwise) and A 2 (clockwise) with the radii satisfying 
0 < Ra 2 < nrniyX) fly and Rax > maxy>ofly, as in 
Fig. a we can apply the series (l 2 | ) and ( fl3] ) along the 
integration path. The odd-fc sum rules are then given as: 

E [ u} k S(F,ui)duj (k > 0,odd), 

mEM 1 /E - E) 

- / oj k S(F,uj)duj (k < 0,odd). 

2 tt* Ja 2 


To evaluate even-fc sum rules, we need to connect A\ 
and A -2 to enclose the positive-energy poles. To this end, 
we consider contour D of Fig. [l] composed of semicircles 
A 1: A 2 connected by straight segments I\ and I 2 on the 
imaginary axis. 
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FIG. 1 . (Color online) The contour D (oriented counter¬ 
clockwise) in the complex-cu plane used to evaluate sum rules. 
The contour consists of two semicircles A\ and A2 (of radii 
Ra 1 and Ra 2 , respectively) and two segments I\ and I2 on 
the imaginary axis. The positive QRPA poles are all located 
between Ra 2 an d Rai- 


In summary, regardless of the moment k , the sum rule 
is given by the integration along D: 


m k (F) 


1 

27ri 


oj k S(F,u)cLj = y £n$ |H^|0)| 2 . 

Jd is>0 


(15) 


However, for certain moments k, some parts of path D 
do not contribute to the sum rule. For odd values of k, 
the contributions from I\ and I 2 cancel each other. Fur¬ 
thermore, for negative fc, application of Jordan’s lemma, 
together with a limit of Ra, —> 00 , allows for the removal 
of the contribution from A\. For positive fc, there is no 
pole at w = 0 , and the limit Ra 2 —> 0 can be taken. 
Table |T] lists the portions of the contour D required for 
each k. Furthermore, for even fc, the contributions from 
I\ and I 2 are identical. Similar contours were considered 
in Ref. [52] to compute energy-weighted sum rules. 


TABLE I. Portions of the contour D required for comput¬ 
ing various sum rules mk . For sum rules with even k, the 
contributions from 7 i and I2 are identical. 


k 

Required portions of D 

negative, even 

A2, Ii, I2 ( Ra 1 — t 00) 

negative, odd 

A2 {Rai —00) 

0 

A 1, A2, 1 1, I2 

positive, odd 

Ai ( Ra 2 0 ) 

positive, even 

Aij /1, I2 ( Ra 2 0 ) 


IV. RESULTS 

A. Numerical checks and benchmarking against 
MQRPA 

To check the FAM approach to sum rules, following 
Refs. J3l [51] we consider the oblate configuration of 
24 Mg computed with the SLy4 [5jS] Skyrme EDF. The 
HFB calculations were carried out with the DFT solver 
HFBTHO 51] in a model space of five harmonic os¬ 
cillator shells by employing a volume pairing with the 
strength of V 0 = —125.20 MeV fm 3 and a 60MeV quasi¬ 
particle energy cutoff. The resulting oblate minimum of 
24 Mg has nonzero pairing in protons and neutrons. The 
small single-particle model space employed makes it pos¬ 
sible to benchmark FAM results against the matrix for¬ 
mulation of QRPA (MQRPA) [55] without any further 
truncation. To compute spatial integrals we used Gauss- 
Hermite (Nq h = 30), Gauss-Laguerre (Nq l = 30), and 
Gauss-Legendre (7V Leg = 30) quadratures. The finite- 
amplitude expansion parameter 77 was set to be 10 -7 , 
and the convergence criterion of the FAM was set such 
that the change of the individual FAM amplitudes from 
the previous iteration should be less than 10 -5 . The inte¬ 
gration along semicircles A\ and A 2 was discretized with 
Na-l and Na 2 points, respectively. In addition, the in¬ 
tegration along I\ was discretized with Nj ± points and 
evaluated using the composite Simpson’s rule. As for 
negative-fc moments, the composite Simpson’s rule was 
applied to the variable 1/oj to describe the divergent be¬ 
havior of integrand around lu = 0. In this particular 
test case, the smallest and largest energy MQRPA poles 
appear at 1.3 MeV and 128.7 MeV, respectively. Conse¬ 
quently, the contour radii were set to Ra 2 = 1 MeV and 
Ra 1 = 200 MeV. In order to systematically assess our 
numerical procedure for different moments k, we used 
the same contour D for all cases, without simplifications 
listed in Table [U 

As far as the external field F is concerned, we consid¬ 
ered the isoscalar (IS) and isovector (IV) monopole (M) 
and quadrupole (Q) operators: 


i= 1 

(16) 

A 


F lYM =Y e e s(T 3i )r^T 3i , 

i=l 

(17) 

piSQ r?Y 2O (0i, ipi), 

i=l 

(18) 

A 


pw Q =^e e ff(r 3i )r 2 l2o(^,^j)T3i, 

(19) 


i=1 


where e e ff(n) = eZ/A and e e s(p) = —eN/A. 

It is worth noting that neutron and proton pairing- 
rotational spurious modes, associated with the break¬ 
ing of the particle number symmetry, are present in 
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the K v = 0+ sector. Fortunately, these modes - gen¬ 
erated by the neutron and proton particle number op¬ 
erators - cannot be excited by particle-hole operators 
(16)-(19). Therefore, the presence of pairing-rotational 


spurious modes does not cause any additional difficulties 
[311143]. 

To begin with, we checked the convergence of the inte¬ 
gral (151 along Ai, A 2 , and I\ with respect of the number 


of integration points. The results are presented in Ta¬ 
bles m for the isoscalar monopole operator. As seen 
in Table |II[ the integrals along Ai are small for negative 
k. Analytically, these values should be zero for negative- 
odd values of k ; hence, nonzero values in Table |H| reflect 
the numerical error of calculations. As far as the posi¬ 
tive k moments are concerned, the convergence is faster 
for odd-fc sum rules. In particular, the convergence for 
k = 1 is excellent, since a 6 -digit accuracy is achieved al¬ 
ready with Na x = 5. The integration along Ai captures 
the total mi and m 3 sum rules; the result in Table [II] 
indicates that these sum rules can be computed very effi¬ 
ciently. Moreover, since the semicircle A\ is located very 
far from the QRPA poles, FAM calculations along A\ 
converge very quickly, typically after 6 iterations. Fur¬ 
thermore, each FAM calculation at a given w along the 
contour is easily parallelizable; this could significantly re¬ 
duce the total computational time, although not so many 
points are required for the convergence of m\ and m 3 . 


Table III shows the convergence of the integral (15) 


along A 2 - This portion of the contour is required for the 
sum rules with negative k. Of most practical importance 
is the inverse energy-weighted sum rule m_i. The value 
of m_ 1 converges here with Na 2 = 16 points. In gen¬ 
eral, as compared to integration along Ai, more FAM 
iterations are required to achieve reasonable convergence 
along A 2 . In the case considered, typically 50 FAM it¬ 
erations are necessary for each ui. When choosing Ra 2 
one has to keep in mind that its value should be smaller 
than the lowest QRPA pole, whose energy is not a pri¬ 
ori known. At the same time, the convergence of FAM 
calculations for negative-A moments deteriorates rapidly 
when Ra 2 gets too close to zero. 

Table [TV] illustrates the convergence along the segment 
Ii on the imaginary axis. As discussed, this integration 
should be nonzero only for even-A moments. The con¬ 
vergence for k = 0 is reached rather slowly, especially 
when compared with the k = 4 and —4 cases. This is 
because the Simpson’s formula used approximates the in¬ 
tegrand with quadratic functions, which is a poor ansatz 
for k = 0 . 

To benchmark our FAM approach, in Table [V] we dis¬ 
play the values of sum rules for the isoscalar and isovec¬ 
tor monopole operators; they are compared with the 
MQRPA results based on the direct evaluation of the 
r.h.s. of Eq. (15). Overall, there is an excellent agree¬ 
ment between the two sets of calculations. This result 
indicates that the proposed FAM technique can be used 
to predict sum rules of interest in model spaces that are 
too large to be treated with MQRPA. The convergence 


of integration along A 2 is not sufficient in the case of 
k = —4; this sum rule is, however, less important than 
other moments discussed. 


B. Thouless theorem for energy-weighted sum rule 

The Thouless theorem [35] gives the relation between 
the energy-weighted sum rule mi(F) for isoscalar F = 

Eti f(fi) or isovector F = fi^hi one-body op¬ 
erators and the expectation value of the double commu¬ 
tator at the ground state 


mi{F) =\([F, [H, .F]]) = ^(1 + k){[F, [T, F]]> 
=( 1 + k)^- j |V f(r)\ 2 p{r)dr , 


( 20 ) 


where T is the kinetic energy operator and n is the en¬ 
hancement factor, which is present in the case when F 
is an isovector operator. The explicit expressions for the 
r.h.s. of Eq. ( |20| ) for the operators (16)-(19) are given in 
Appendix [A] 

The theorem is exact when both the ground state and 
the excited states are many-body shell-model states, and 
has been proven for HF+RPA [35], HFB+QRPA [36] . 
and second RPA m In the case of HF+RPA, expression 
(20 1 also holds for the Skyrrne force due to the ^-character 
of the momentum-dependent terms [231 S31 . However, 
as pointed out in Ref. [33] . the theorem has not been 
proven for a generalized EDF, which is not explicitly re¬ 
lated to an effective interaction. Deviations from relation 


( 20 ) can be caused by, e.g., different assumptions about 


particle-hole and pairing channels, the Slater approxima¬ 
tion to the Coulomb exchange term, approximations to 
spin-orbit and tensor terms [53] . and generalized density 
dependence [59H6T] . To the best of our knowledge, the 
Thouless theorem has not been proven in the case of gen¬ 
eralized EDFs. 


In the following, we refer to the value (20) as the “HFB 


value” of the energy-weighted sum rule. In Table VI the 
energy-weighted sum rules obtained in HFB and FAM 
are compared for different model spaces given by N s h. In 
a small model space of N s h = 5, the difference between 
FAM and HFB values is non-negligible but quickly be¬ 
comes small with N^. This can be attributed to a poor 
representation of the operator F in small basis spaces, re¬ 
sulting in an error on the derivative of the function f(f) 
in Eq. (20). In spite of the fact that the SLy4 EDF com¬ 


bined with volume pairing cannot be related to a force, 
the numerical test in Table El demonstrates that the 
Thouless theorem provides a reasonably good approxi¬ 
mation to the value of the sum rule mi for the Skyrrne 
EDF. 


In the notation of Ref. 


the time-odd part of the 
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TABLE II. The real part of the integral (151 for —4 < k < 4 and F lb (in Me v e“ fm 4 ) along semicircle Ai with Ra 1 = 
200 MeV. The integral was discretized with Na 1 points. The numbers in parentheses denote powers of 10. 


N Al 

k = —4 

CO 

1 

II 

-56 

k = -2 

k = -1 

k = 0 

k = 1 

k = 2 

CO 

II 

-*6 

k = 4 

2 

-9.0(-9) 

-2.6(-6) 

-3.8(-4) 

-2.7(-3) 

14.4510 

4195.31 

608575 

4318446 

-23121542422 

3 

9.0(-9) 

6.8(-8) 

-1.7(-4) 

l.l(-4) 

13.8328 

4199.78 

574147 

4338677 

-10598058261 

4 

3.3(-9) 

-2.6(-9) 

-1.4(-4) 

-5.7(-6) 

13.5751 

4199.39 

562629 

4336463 

-8502934038 

5 

2.5(-9) 

2.2(-10) 

-1.3(-4) 

4-1 (-6) 

13.4599 

4199.44 

557421 

4336545 

-7722808528 

7 

2.0(-9) 

7.1(-12) 

-1.2(-4) 

2.9(-6) 

13.3593 

4199.44 

552937 

4336634 

-7119788306 

9 

1.9(-9) 

2.1 (-12) 

-l.l(-4) 

2.9(-6) 

13.3181 

4199.44 

551106 

4336643 

-6889795483 

10 

1.8(-9) 

1.9(-12) 

-l.l(-4) 

2.9(-6) 

13.3061 

4199.44 

550575 

4336644 

-6824728295 

11 

1.8(-9) 

1.9(-12) 

-l.l(-4) 

2.9(-6) 

13.2972 

4199.44 

550183 

4336638 

-6777106568 

12 

1.8(-9) 

2.0(-12) 

-l.l(-4) 

2.9(-6) 

13.2905 

4199.44 

549885 

4336648 

-6741178415 

101 

1.6(-9) 

1.9(-12) 

-l.l(-4) 

2.9(-6) 

13.2556 

4199.44 

548341 

4336643 

-6558793102 

201 

1.6(-9) 

1.9(-12) 

-l.l(-4) 

2.9(-6) 

13.2552 

4199.44 

548325 

4336643 

-6556876296 

301 

1.6(-9) 

2.0(-12) 

-l.l(-4) 

2.9(-6) 

13.2551 

4199.44 

548322 

4336644 

-6556517206 


TABLE III. Similar to Table [TT] but for the integration along A 2 for several values of Na 2 with Ra 2 = 1.0 MeV. 

Na 2 

k = —4 

k = - 3 

k = - 2 

k = -1 

O 

II 

-56 

k = 1 

k = 2 

CO 

II 

k = 4 

5 

-1.20929 

0.0372710 

3.26165 

5.00438 

3.23108 

1.0(-3) 

-1.22915 

2.2(-3) 

0.996926 

11 

-1.06784 

0.0369978 

3.21892 

5.00387 

3.18902 

2.7(-5) 

-1.09039 

4.1(-5) 

0.691152 

15 

-1.05236 

0.0369925 

3.21386 

5.00386 

3.18407 

7.1(-6) 

-1.07507 

3.9(-6) 

0.663895 

16 

-1.05021 

0.0369910 

3.21315 

5.00385 

3.18338 

4.1(-6) 

-1.07293 

-7.8(-7) 

0.660178 

21 

-1.04368 

0.0369929 

3.21099 

5.00385 

3.18127 

6.2(-6) 

-1.06646 

-6.4(-7) 

0.649056 

31 

-1.03883 

0.0369918 

3.20937 

5.00385 

3.17968 

6.1(-6) 

-1.06164 

9.9(-7) 

0.640898 

51 

-1.03625 

0.0369916 

3.20851 

5.00385 

3.17884 

5.6(-6) 

-1.05908 

1.0(-6) 

0.636594 

101 

-1.03512 

0.0369918 

3.20813 

5.00385 

3.17847 

5.9(-6) 

-1.05797 

7.2(-7) 

0.634728 


TABLE IV. Similar to Table [TT] but for the integration along I\ for several values of Ni ± with Ra 1 = 200 MeV and Ra 2 = 
1.0 MeV. 


Ni, 

1 

II 

-Sfi 

CO 

1 

II 

-56 

k = -2 

k = -1 

O 

II 

-56 

k = 1 

k = 2 

CO 

II 

-56 

k = 4 

10 

0.523854 

-6.6(-8) 

-1.472654 

1.3(-7) 

61.4767 

-2.0(-6) 

-208884 

2.3(-2) 

3356473886 

30 

0.523849 

1.2(-6) 

-1.478626 

-1.6(-6) 

61.7093 

-6.0(-7) 

-208523 

-4.2(-3) 

3356786248 

50 

0.523850 

-l.l(-6) 

-1.477671 

1.4(-6) 

61.7024 

-2.2(-6) 

-208522 

3.3(-2) 

3356787504 

100 

0.523850 

-1.5(-6) 

-1.477427 

2.0(-6) 

61.7044 

-1.4(-6) 

-208522 

2.8(-2) 

3356787251 

200 

0.523850 

-8.7(-7) 

-1.477417 

1.2(-6) 

61.7052 

-1.8(-6) 

-208522 

2.2(-2) 

3356787739 


TABLE V. Sum rules (in MeV fe e 2 fm 4 ) for the isoscalar and isovector monopole operators calculated with MQRPA and FAM. 
The FAM calculations were performed by using Na 1 = 301, Na 2 = 101, and Ni 1 = 200 integration points. 



k = —4 

k = - 3 

k = - 2 

k = -1 

O 

II 

-56 

k = 1 

k = 2 

k = 3 

k = 4 

MQRPA(ISM) 

0.013077 

0.037185 

0.253118 

5.00072 

139.825 

4200.82 

131368 

4342358 

157906069 

FAM(ISM) 

0.012579 

0.036992 

0.253186 

5.00385 

139.844 

4199.44 

131277 

4336644 

157058272 

MQRPA(IVM) 

0.00063616 

0.00273872 

0.07120949 

2.78540 

113.908 

4735.03 

199525 

8527358 

370625216 

FAM(IVM) 

0.00043157 

0.00275227 

0.07133510 

2.78615 

113.908 

4734.30 

199510 

8524830 

368643941 
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TABLE VI. The energy weighted K n = 0 + sum rule (in MeV e 2 fm 4 ) for the operators (161-1191 at the oblate minimum of 
24 Mg as a fu nct ion of N s h- The FAM values were obtained by taking Ra-i = 200 MeV and Nai = 12; they are compared to 
HFB values (201. The results without time-odd terms except for the current-current coupling (C 3 t ^ 0 and Ct(po) = Cl a = 
Cj* = Cj = 0) (a), and with the full time-odd functional except for the current-current and kinetic spin-spin couplings 
(C{ = Cl = 0, C a (po) ^ 0, Cl s ^ 0, and C t v3 ^ 0) (6), obtained with A s h=20, are also listed. 


N ah 

FAM(ISM) 

HFB(ISM) 

FAM(IVM) 

HFB(IVM) 

FAM(ISQ) 

HFB(ISQ) 

FAM(IVQ) 

HFB(IVQ) 

5 

4199.44 

4303.67 

4734.25 

4752.37 

762.638 

767.933 

848.110 

845.235 

10 

4524.39 

4502.75 

4970.34 

4940.08 

779.019 

775.724 

852.995 

849.015 

15 

4521.39 

4523.80 

4958.02 

4960.52 

776.587 

776.161 

849.482 

849.116 

20 

4530.01 

4529.46 

4966.98 

4966.01 

777.425 

776.832 

850.145 

849.747 

20 (a) 

4530.07 

- 

4966.98 

- 

777.506 

- 

850.132 

- 

20 (6) 

5297.64 

- 

5298.46 

- 

905.461 

- 

905.441 

- 


Skyrme EDF reads: 


£odd = E [ C tiP°) s2 t + C ^ s t ' Aa t 

t=0,l 

+ C(j 2 + C? j s t ■(Vxj,) + Cl s t ■ T t \ . 


( 21 ) 


By taking the Skyrme interaction as a starting point, 
the time-odd and time-even coupling constants of the 
Skyrme EDF are related to each other. That is, by fix¬ 
ing time-even coupling constants, the time-odd part be¬ 
comes also determined. This choice also guarantees the 
EDF’s gauge invariance |63| . In the EDF picture, how¬ 
ever, the time-odd coupling constants could be treated 
as independent parameters, where some of them can be 
constrained by local gauge invariance [B2] [M] . With local 
gauge invariance assumed and tensor terms excluded, the 


last term of Eq. (21), proportional to (7/ , vanishes. In 


standard HFB calculations for even-even nuclei, the time- 
odd fields do not contribute because of time-reversal sym¬ 
metry; hence, the time-odd part (21) does not affect the 


HFB value (201. However, when time-reversal symme¬ 


try becomes broken, as in the case of FAM calculations, 
time-odd terms become active. 


As shown in Table VI the inclusion of the current- 
current term Cjj 2 is necessary in the FAM to recover 
the HFB value of the energy-weighted sum rule of the 
monopole and quadrupole operators. This indicates that 
the gauge invariance of the term pr — j 2 should not be 
broken when applying the Thouless theorem to QRPA 
sum rules. Other terms in the time-odd functional do 
not impact the energy-weighed sum rule. Local gauge 
invariance also couples the C t Vj and terms, but the 
numerical results demonstrate that these do not affect 
the energy-weighted sum rule. 


C. Dielectric theorem for the inverse 
energy-weighted sum rule 

The dielectric theorem connects the inverse-energy- 
weighted sum rule (related to nuclear polarizability) with 
the constrained potential energy surface. This theorem 
was proposed in Refs. [20) [23] for the HF case, and has 


been proven in the HFB framework in Ref. m- Based on 
this theorem, the inverse energy-weighted sum rule m_i 
can be obtained from the curvature of the total energy £ 
at equilibrium: 


m -' {p) = \ lv (A) 


A=0 

1 d (0(A) | F |0(A)) 


d\ 


( 22 ) 


A=0 


where the constrained HFB state |</>(A)) is obtained by 
minimizing the total Routhian con tain ing a linear con¬ 
straint —A F. We use the relation (22) to compute the 


m_i sum rule. The derivative is evaluated with a finite 
difference of A A = 0.0001 MeV e -1 fm -2 . The resulting 
rri-i values are compared with those from the FAM in 


Table |VII| A good agreement is found already in a small 
model space (N s h = 5) where to_i is not fully converged, 
indicating that the dielectric theorem works well, inde¬ 
pendently of the size of the model space. This finding is 
consistent with the proof of Ref. m , which applies to an 
arbitrary size of quasiparticle space. 


D. Example of systematic calculations 

As an illustrative example, we discuss the energy- 
weighted A' 71 ’ = 0 + sum rules in the shape transitional re¬ 
gion of 142_152 Nd and 144_154 Sm. The calculations were 
carried out by using the SLy4 EDF parameterization with 
a volume pairing strength V n = V p = —190 MeV fm 3 in 
the model space of N s h = 20 oscillator shells. The pair¬ 
ing strength was adjusted to reproduce the experimental 
proton pairing gap of 1.23 MeV in 142 Nd. In this realis¬ 
tic calculation we use Nq h = Agl = 40 and AL eg = 80, 
which are the recommended values based on recent anal¬ 
ysis [54] . The FAM contour integration was carried out 
using a semicircle with Ra 1 = 200 MeV, discretized with 
Nai = 12 points. 

Table IVIIII summarizes the results. The calculated 
ground state properties show a gradual spherical-to- 
deformed shape transition with increasing neutron num¬ 
ber. Moreover, in some of the isotopes we predict pairing 























TABLE VII. Inverse energy-weighted sum rule (in MeV 1 e 2 fm 4 ) computed using the dielectric theorem (HFB) and the FAM 
for various sizes of the model space given by N s \,. FAM calculations were performed using Na 2 = 22 and Ra 2 = 1.0 MeV. 


N ah 

FAM(ISM) 

HFB(ISM) 

FAM(IVM) 

HFB(IVM) 

FAM(ISQ) 

HFB(ISQ) 

FAM(IVQ) 

HFB(IVQ) 

5 

5.00385 

5.00375 

2.78615 

2.78614 

4.44830 

4.44765 

0.798680 

0.798680 

10 

11.2033 

11.2102 

5.09467 

5.09671 

5.21547 

5.21586 

1.07516 

1.07524 

15 

12.4930 

12.5009 

5.71677 

5.71960 

5.31250 

5.31268 

1.12916 

1.12910 

20 

12.9506 

12.9634 

6.06842 

6.07304 

5.35499 

5.35730 

1.15744 

1.15771 


TABLE VIII. Isoscalar monopole and quadrupole energy-weighted K w = 0 + sum rules in units of MeV e 2 fm 4 computed with 
the FAM and the HFB techniques for 142_152 Nd and 144 ~ 154 Sm isotopes. The quadrupole deformation /3, neutron and proton 
pairing gaps (A n and A p , respectively), and total rms radius \J{r 2 ) are also shown. 


0 

A n (MeV) 

Ap(MeV) 

V / < r ^>( fm ) 

HFB(ISM) 

FAM(ISM) 

HFB(ISQ) 

FAM(ISQ) 

142 Nd 

0.0 

0.00 

1.21 

4.92 

50497 

50724 

10046 

10068 

144 Nd 

0.09 

0.49 

1.09 

4.95 

50453 

50647 

10606 

10626 

146 Nd 

0.15 

0.55 

1.00 

4.99 

50590 

50769 

11042 

11062 

148 Nd 

0.21 

0.00 

0.93 

5.03 

50788 

50936 

11412 

11429 

150 Nd 

0.31 

0.64 

0.00 

5.11 

51667 

51806 

12287 

12301 

152 Nd 

0.32 

0.00 

0.00 

5.14 

51649 

51762 

12375 

12383 

144 Sm 

0.0 

0.00 

1.10 

4.94 

53635 

53873 

10670 

10693 

146 Sm 

0.06 

0.55 

1.08 

4.97 

53492 

53712 

11048 

11069 

148 Sm 

0.16 

0.56 

1.07 

5.01 

53754 

53957 

11770 

11792 

150 Sm 

0.21 

0.16 

0.93 

5.06 

53979 

54145 

12189 

12207 

152 Sm 

0.28 

0.57 

0.69 

5.11 

54474 

54646 

12768 

12786 

154 Sm 

0.32 

0.09 

0.65 

5.16 

54707 

54849 

13071 

13084 


collapse. For that reason, the chosen set of nuclei is rep¬ 
resentative of a realistic situation encountered in global 
surveys across the nuclear landscape, where deformations 
and pairing may vary rapidly as a function of proton and 
neutron number. 

The energy-weighted sum rules computed with the 
FAM agree well with the HFB expressions of Appendix[A| 
This agreement holds regardless of nuclear shape or pair¬ 
ing. As expected, the energy-weighted sum rule for the 
isoscalar monopole operator increases with N in the re¬ 
gion of the shape transition; this is attributed to an 
increase of the root mean square radius with deforma¬ 
tion. Similarly, the isoscalar quadrupole operator in¬ 
creases even more rapidly with increasing quadrupole de¬ 
formation. 

Next, we consider the energy weighted sum rules in 
constrained HFB states. The constrained HFB poten¬ 
tial energy curve as a function of quadrupole moment 
was obtained using the quadratic constraint. The con¬ 
tribution from the quadratic constraining potential was 
included consistently to the residual field in the FAM. 
This kind of calculation represents local QRPA on top of 
constrained HFB |65| : it contains dynamical information 
about non-equilibrium configurations in the deformation 
space. 

The energy-weighted sum rule of the isoscalar 
quadrupole operator as a function of quadrupole defor¬ 
mation is shown in Fig. [2] The sum rule increases mono- 
tonically with /3 and agrees very well with HFB values. 


This, together with results presented in Table |VIII[ in¬ 
dicates that the Thouless theorem provides a good ap¬ 
proximation to the energy weighted sum rule within the 
Skyrme-EDF picture, which is not based on the underly¬ 
ing Hamiltonian. 



P 

FIG. 2. (Color online) Energy-weighted sum rule of the 
isoscalar quadrupole operator in 142 Nd obtained in the HFB 
(solid line with asterisks) and FAM (dashed line with open 
circles) frameworks as a function of quadrupole deformation 
/? for the constrained HFB solutions. 

In passing, we should note that when departing from 
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the HFB minimum, there is a possibility of imaginary 
energy QRPA solutions; in such cases, a pair of QRPA 
poles would appear on the imaginary axis. Although one 
expects no contribution to odd-fc sum rules from such 
a pair, a careful consideration needs to be given to the 
choice of integration contour in the FAM. A general ex¬ 
tension of the complex FAM formalism to the case of local 
QRPA will be an interesting avenue for future studies. 


V. CONCLUSIONS 

We propose an efficient formalism to compute sum 
rules by using the contour integration formalism within 
the complex-energy finite-amplitude method. In partic¬ 
ular when the order of the moment is odd, the obtained 
expression becomes extremely simple, as the sum rules 
appear as expansion coefficients of the Laurent series of 
the response function. The new formalism has been suc¬ 
cessfully benchmarked against the matrix diagonalization 
method of QRPA. 

We compare the energy-weighted sum rule obtained in 
the FAM with those based on the Thouless theorem. Al¬ 
though the double commutator cannot be evaluated for 
general energy density functionals that are not based on a 
Hamiltonian, the numerical results indicate that the the¬ 
orem provides a very good approximation to m 1 when a 
large model space is employed and local gauge symmetry 
of the EDF is satisfied. The inverse-energy-weighted sum 
rule was compared with the constrained HFB result us¬ 
ing the dielectric theorem, and a perfect agreement was 
obtained regardless of the model space. 

Our results suggest that sum rules can be computed 
efficiently in the FAM even in cases when other methods 
are not easily available (e.g., the Thouless theorem can¬ 
not be applied or constrained calculations cannot be car¬ 
ried out because of self-consistent symmetries assumed). 
The extension of the formalism to non-Hermitian opera¬ 
tors is also straightforward, as it has already been applied 
to the beta-decay rates with pn-FAM [44] . 

The FAM approach to sum rules promises to add 
new functionality to the EDF optimization framework 
of Refs. [SHH] as it will allow adding new kinds of data 
on multipole- and charge-exchange strength to the set 
of fit-observables defining the objective function. The 
new FAM technique can be very useful when studying 
the nuclear response to non-trivial operators such as the 
nuclear Schiff moment, which is closely related to the 
isoscalar dipole operator mm- 
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Appendix A: Thouless theorem for monopole and 
quadrupole operators 


According to the Thouless theorem (201, the energy- 
weighted sum rule for isoscalar monopole and quadrupole 
operators of an axially-deformed nucleus are: 


mi (ISM) 
toi(ISQ) 



2 h 2 

—Mr 2 ) 

2 m X ' 


h 2 5 

2TO 2-7T 


Mr 2 ) 



(Al) 

(A2) 


where (r 2 ) is the total rms squared radius and /3 is the 
mass quadrupole deformation parameter: 


P / (3 ~ 2 ~ r ^p^ dr - ( A3 ) 

For isovector operators, there appears an enhancement 
factor 


K = ^{c%-ciy 


Mf(r)\ 2 p n (r)p p (r)dr 
I Mf(r)\ 2 p(r)dr 


(A4) 


where - in the notation of Ref. [62] - CJ is the coupling 
constant of the term p t r t in the EDF. The expressions for 
the isovector monopole and quadrupole operators are: 


mi(IVM) ==4e2 ^ n ^f [ z ( r 2 )n + N ( r 2 ) p ] (1 + kivm), 

(A5) 

«ivm = 7 (t( C 'o - CT)^py J r 2 p n (r)p p {r)dr , 

(A6) 


and 
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mi(IVQ) 


2 h 2 NZ 5 
2m A 2 2 tt 




where subscripts n/p indicate neutron/proton expecta¬ 
tion values respectively. 


+ N(r 2 } p 



(1 + Kivq), 


(A7) 


«IVQ 


8m 

H 2 


(C T 0 


ci) - 

2 A(r 2 ) 



x J (3z 2 


r 2 )pn(r)p p (r)dr, 


(A8) 
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